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Some  Probability  Models  Related 
To  Classical  Occupancy  With  Possible 
Applications  To  Tactical  Analysis 

by 
Bruno  0.  Shubert 


Abstract 
In  this  report  three  types  of  probability  distributions  associated 
with  random  search-type  models  are  derived.  The  distributions  are  closely 
related  to  the  classical  occupancy  problem.  The  first  two  types  are  uni- 
variate and  are  applicable  to  simple  random  search  models  such  as  mine 
hunting  or  one-kind  mine  sweeping.  The  third  type  is  multivariate  allowing 
for  various  mine  ship  count  settings. 


1 .   Introduction 

In  this  report  we  investigate  some  discrete  probability  models  which 
may  be  useful  in  search  and  related  situations.  The  models  are  closely 
related  to  the  classical  occupancy  problem  of  probability  theory  and  may  in 
fact  be  viewed  as  generalizations  of  it. 

The  models  are  discrete  in  the  sense  that  there  are  m  discrete 
locations  -  referred  to  as  eel  Is  -  possibly  containing  the  objects  of  the 
search  -  referred  to  as  bal Is.  The  search  itself  is  random  in  the  sense 
that  the  cell  to  be  searched  is  selected  at  random,  i.e.  with  equal 
probability  among  all  the  cells  involved. 

Two  univariate  models  and  a  multivariate  extension  of  the  first  of  the 
two  are  considered.  The  detailed  description  of  the  first  univariate 
model  -  called  here  Random  Search-Variant  1  -  is  as  follows: 

There  are  m  cells,  y   of  them  initially  empty  and  y,  of  them 
containing  a  single  ball  each.  The  search  consist  of  repeated  independent 
trials,  where,  at  each  trial,  a  cell  is  selected  with  probability  1/m  .  The 
selected  cell  is  examined,  and  if  it  contains  a  ball  the  ball  is  found  with 
probability  p  .   If  a  ball  is  found,  it  is  removed  from  the  cell  so  that 
the  cell  is  now  empty.  Since  the  trials  leading  to  a  cell  selection  are 
assumed  independent  the  same  cell  can  be  selected  again  regardless  of  the 
results  of  previous  trials.  Thus  a  cell  can  be  selected  even  if  a  ball  was 
found  and  removed  from  it  at  some  previous  trial. 

This  model  is  investigated  in  Section  3.  We  find  the  probability  dis- 
tribution of  the  number  of  balls  found  in  n  trials,  the  number  of  trials 
needed  to  find  the  k-th  ball  as  well  as  moments  of  these  distributions. 

As  mentioned  earlier,  this  model  can  be  rephrased  as  an  occupancy 
problem  of  distributing  n  balls  into  m  cells.   In  fact,  if  y  =  0  and 


p  =  1  it  is  identical  with  the  classical  occupancy  problem  -  merely  iden- 
tify the  number  of  balls  found  in  our  model  with  the  final  number  of  nonempty 
cells  in  the  occupancy  problem. 

The  second  univariate  model  -  Random  Search  Variant  2  -  differs  from 
Variant  1  only  in  the  rule  for  selecting  cells  at  each  trial.  This  time 
we  assume  that  the  cell  in  which  a  ball  has  been  found  and  removed  is  no 
longer  a  candidate  for  further  selections.  Thus,  the  number  of  cells  to  be 
searched  is  reduced  by  one  each  time  a  ball  is  found.  The  corresponding 
distributions  of  the  number  of  balls  found,  the  number  of  trials  needed  to 
find  the  k-th  ball  and  their  moments  are  described  in  Section  4.  It  may  be 
interesting  to  notice  that,  with  p  =  1,  this  model  is  a  hybrid  of  a  classi- 
cal sampling  with  and  without  replacement.  Think  of  the  m  cells  as  m 
balls,  y   of  them  white  and  y,  black.  Then  draw  n  balls  replacing 
the  white  ones  but  not  replacing  the  black  ones,  and  ask  about  the 
distribution  of  the  number  of  black  balls  drawn. 

The  last  section  of  this  report  deals  with  a  multivariate  extension 
of  the  Variant  1  model.  The  extension  consists  in  assuming  that  there  are 
now  s  different  types  of  balls  involved  and  each  time  a  ball  of  type 
j  >  0  is  found  it  is  replaced  by  a  ball  of  type  j  -  1  or  removed  if 
j  -  1  =  0  .  The  model  parameters  are  now  nonnegative  integers  y  , ..., 
y  ,  where  y  +...+  y  =  m  ,  y   is  the  number  of  empty  cells  and  y.  is 
the  number  of  cells  with  balls  of  type  j  each  initially.  For  simplicity, 
only  the  value  p  =  1  is  considered  here,  i.e.  if  a  cell  containing  a  ball 
is  selected  the  ball  is  found  (and  its  type  reduced)  with  probability  one. 

The  cell  selection  process  is  the  same  as  in  Variant  1.  That  is,  cells 
are  selected  at  random  among  the  m  cells  in  independent  trials  so  that  the 
same  cell  can  be  selected  repeatedly  regardless  of  the  ball  type  it  contains. 


Again,  this  model  can  be  rephrased  as  an  occupancy  problem  of  distributing 
n  balls  into  m  cells  and  interpreting  a  cell  with  ball  type  j  as  a  cell 
containing  s  -  j  balls  for  j  >  0  or  a  cell  containing  at  least  s  balls 
for  j  =  0  .  However,  classical  occupancy  problems  invariably  assume  that 
that  initially  all  cells  are  empty,  i.e.  that  y  =  m  in  our  notation. 
Therefore,  classical  results  are  not  directly  applicable  to  our  model. 

In  Section  5  we  find  the  first  two  (joint)  moments  of  the  resulting 
numbers  of  balls  type  j  =  Q,...,s  after  n  trials  as  well  as  of  the  num- 
bers of  each  ball  type  found.  An  expression  for  the  marginal  distribution 
of  the  number  of  balls  type  j  is  also  found  -  although  the  complexity  of 
this  expression  leaves  some  doubts  about  its  possible  uses. 

In  this  report,  no  attempt  is  made  to  study  various  asymptotic  distri- 
bution resulting  in  letting  the  parameters  of  the  models  increase  to  infin- 
ity in  some  way.  Although  this  may  yield  a  considerable  simplification,  for 
instance  by  a  Poisson  distribution  in  the  case  of  Variant  1,  in  most  of  the 
applications  intended  here  the  values  of  the  parameters  are  typically  small. 
Still,  asymptotic  results  may  perhaps  warrant  future  investigation. 

Likewise,  no  algorithms  to  actually  compute  numerical  values  for  the 
derived  quantities  are  presented  in  this  report  -  mainly  because  of  a  limited 
time  and  each  of  adequate  computational  facilities  during  preparation  of 
this  report.  This,  again,  is  left  for  a  possible  future  work. 


2.  Possible  Applications 

Although  the  purpose  of  this  report  is  to  develop  model  building  tools 
rather  then  real  world  models  themselves,  a  few  possibilities  will  be 
mentioned  briefly. 

A  rather  straightforward  application  of  the  two  univariate  cases  would 
apply  when  one  indeed  has  m  distinct  geographical  locations  out  of  which 
y,  contain  at  most  one  target  each.  Searching  these  locations  at  random 
would  imply  an  absence  of  any  particular  search  strategy  or  any  further  in- 
formation about  the  locations  of  these  targets.  Such  a  model  may,  for  in- 
stance, serve  as  a  basis  against  which  an  efficiency  of  various  search 
strategies  can  be  compared.  In  the  case  when  the  targets  are  distributed 
in  a  continuous  geographical  region,  the  discrete  cells  may  be  defined  arti- 
ficially by  partitioning  the  region.  The  actual  size  of  a  cell  must  be 
chosen  to  satisfy  the  conflicting  requirements  of  the  size  of  the  searcher's 
detection  region  and  the  assumption  of  a  single  target  in  each  cell.  This, 
of  course,  is  the  problem  encountered  any  time  a  continuous  model  is  being 
discretized.  The  choice  between  Variants  1  and  2  is  then  dictated  by  the 
nature  of  the  targets  and/or  searcher.  For  example,  if  the  targets  are  sta- 
tionary, Variant  2  is  more  appropriate.  Typical  example  here  could  be 
underwater  mine  hunting.   If  the  targets  are  mobile  and  can  randomly  re- 
arrange between  each  search  trial,  Variant  1  could  be  used.  Another  possi- 
bility leading  to  Variant  1  is  when  the  searcher  (in  the  general  sense  of 
this  term)  has  no  control  over  the  cell  selection.  Typical  application  in 
this  case  is  of  course  an  artillery  coverage  model  -  delivering  n  shells 
to  a  region  partitioned  into  m  cells  out  of  which  y,  contain  a  target. 
Similar  applications,  not  involving  an  actual  physical  search,  are  con- 
ceivable. For  instance,  in  reliability  one  can  have  an  equipment  with  two 


kinds  of  components  -  one  with  potentially  unlimited  supply  of  spares  and 
the  other  with  any  a  small  number  y,  of  spares. 

The  possible  applications  of  the  multivariate  version  are  even  more 
versatile  depending  on  the  interpretation  of  the  ball  type.  Consider,  for 
instance,  a  mine  sweeping  model  with  a  ball  type  being  the  ship  count  the 
mine  (ball)  is  initially  set  at.  Alternately,  consider  a  mine  field  with 
cells  corresponding  to  paths  through  the  mine  field  and  ball  type  corre- 
sponding to  the  number  of  mines  along  a  path.  The  cells  would  be  separated 
from  each  other  by  the  actuation  width  of  a  mine  in  this  application.  Yet 
another  application  is  again  the  artillery  coverage  with  multiple  hits 
required  to  destroy  a  target. 

Our  final  remark  concerns  the  choice  of  the  model  parameters.  They 
can,  of  course,  be  left  as  true  parameters  to  investigate  their  influence 
on  the  quantity  of  interest.  But  they  themselves  can  be  considered  random 
variables  resulting  in  various  mixtures  of  the  probability  distributions 
involved.  For  instance,  the  number  y-,  of  targets  could  be  a  Binomial  ran- 
dom variable  thus  modeling  the  effect  of  a  random  appearance  of  targets. 
Even  more  realistically,  perhaps,  the  ship  counts  in  the  mine  sweeping  ap- 
plication can  be  considered  set  at  random  among  some  given  range  of  ship 
counts.  Such  assumptions  could  sometimes  even  simplify  the  analysis  as 
shown  e.g.  at  the  end  of  Section  5. 


3.  Random  Search  -  Variant  1 

For  the  sake  of  reference  let  us  restate  the  assumption  of  the  model. 
We  have  m  cells,  y   of  them  empty  and  y,  of  them  containing  a  single 
ball  each.  Cells  are  being  selected  repeatedly  with  equal  probability  1/m 
in  independent  trials.  If  a  selected  cell  contains  a  ball  the  ball  is  re- 
moved with  probability  p  >  0  and  remains  there  with  probability  q  =  1  -  p 
independently  of  the  results  of  previous  trials.  Thus,  in  Variant  1,  the 
same  cell  can  be  selected  again  regardless  of  the  event  that  a  ball  may  have 
been  removed  from  it  before. 

The  parameters  of  the  model  are 

y0  .  y-|  >  and  p  , 
where  y  ,  y-,  are  nonnegative  integers  and  0  <  p  <  1  .  The  letter 
m  =  y0  +  y]  . 


Let  K(n)  be  the  number  of  balls  obtained  at  the  conclusion  of  the 
n-th  trial.  Clearly  0  ^  K(n)  £min{n,y-,}  and 


(3.1) 


K(n) 


with  probability  1  - 


y-,  -  K(n) 


K(n+1)  = 


y,  -  K(n) 
K(n)  +  1    with  probability    — - p 


for  n  =  0,1,...  with  K(o)  =  0 
Hence  calling 


Pn(k)  =  P(K(n)  =  k)  , 


we  have  the  recurrence 


(3.2) 


Wk>  = 


1  - 


yi  -  fc  1 


P  (k) 
nv  ' 


y,  -  k  +  1 

P  Pn(k-1)  , 
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n  =  0,1,...  ;  with  Pn(-1)  =  0  and  P  (0)  =  1  ,  P  (k)  =0  for  k  >  0  . 

This  can  be  used  to  recursively  evaluate  the  distributions  P  (k)  for 
given  values  of  the  parameters  m,  y,  and  p  .  However,  it  is  also 
possible  to  obtain  an  explicit  formula,  for  instance  by  employing  the 
generating  functions 


(3.3) 


Mt)  =  I      t"  P  (k)     ,  k  =  0,1,.  ...y. 

K     n=0     n 


From  (3.2)  with  k  =  0  we  have 


1  -  —  p 


Pn(0)  , 


and  using  pn(0)  =  1 


P  (0) 
nv  ' 


yl 
1  --FP 


,  n  =  0,1, 


Hence 


(3.4) 


GQ(t)  - 
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1  -  t 


,   yl 

1  'IP 


From  (3.2)  with  k  >  0  by  multiplying  by  t    and  summing  over  n  =  0,1, 
we  obtain 


Gk(t)  +  t 


y1  -  k  +  1 


P  Gk_-,  (t) 


whence  solving  for  G. (t)  and  iterating 


,U)M 


(3.5) 
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,  k  =  0,... ,y]  , 


(k)     yr  -1 

where  y|  '  =  -> ,  t ,  .  Since  the  k  +  l  roots  t..  =  a,.  , 


J    J 


yi  -  j 

a.  =  1 —  p  ,  of  the  denominator  are  all  distinct  we  can  use  the 


partial  fraction  expansion 
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yl   "  J 

With     a.  =  1 —  p     the  coefficients  are 
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Upon  substitution  into  (3.5)  we  have 


Gk(t) 


yjfc)  tk   V   (.l)k-J 
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whence  expanding  (1  -  ta.)~  into  a  geometric  series 
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The  coefficients  of  this  expansion  are  the  desired  probabilities,   i.e. 


(3.6) 
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Remark 


k  =  0,1,... ,y1    ;   n  =  0,1  ,.. 


(3.7)  Since  I     (-l)jM(x  +  k  -  j)n  =  Akxn 

j=0  ^Jj 


with     a     the  forward  difference  operator  one  can  also  write 


P  (k)   = 
nv    ' 


'SfP  ^k 


r  >  n 

m 

l?-yi 


With  y,  =  m  and  p  =  1  this  is  the  classical  occupancy  distribution 
([Fl],  p.  58)  although  the  generalizations  with  y,  <  m  or  p  <  1  can 
also  be  found  in  the  literature  ([JK],  p.  124,  140).  Note  that  the  formula 
(3.6)  gives  automatically  P  (k)  =  0    for    k  >  n  . 

Next,  let  us  look  at  the  moments  of  this  distribution.  The  expectation 
u  =  E[K(n)]  (and  higher  moments  as  well)  can  be  derived  from  the  recurrence 
(3.2)  by  conditioning.  For  instance  taking  conditional  expectation  we 
obtain 
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However,  one  can  obtain  all  factorial  moments 


»lr)   -   E[K(n)(K(n)  -  l)...(K(n)  -  r  +  1)]  , 
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directly  from  (3.7)  since 
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and  linearity  of  the  difference  operator.  Equivalently 
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whence  again  the  expectation 
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Waiting  time  distributions: 

Consider  next  the  random  variable  W.  defined  to  be  the  trial  number 
at  which  the  k-th  ball  is  found.  We  have 


Wk=T1+...+Tk, 


k.  =  1 ,...  ,y1  , 


where  T.  is  the  number  of  trials  to  find  the  j-th  ball  counted  from  the 
the  trial  at  which  the  (j-l)-th  ball  was  found.  Clearly, 

y1  -  J  +  i 


T.  is  geometric  with  the  parameter 
independent.  Thus  calling 
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where  G  is  the  generating  function  (3.3). 
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k-l 


1  --TT-P 


>+k-l 


and  thus 


P(W.    =  n)  =  ^]       T     (-1)k-1-J'(k"r 
k  n  [k  J    .LQ         ;  [  j  J 


i   yi : J 
1  -—i—p 


n-l 


=    f*] 


v 

k 


k-l 
I     (-DJ 

J-o 


fk-lHm  ,       , 


n-l 


n  =   k,   k+1 ,...    ;     k  =  1 ,. . .  ,y]    . 


Again,  this  can  be  reexpressed  as 


p<wk  ■  »>  "  fef  "  ? 


k-l  fm  f"1 

»        p-*l 


The  expectation  and  the  variance  can  be  gotten  most  easily  by  using 
the  fact  that     W.      is  a  sum  of  independent  geometric  random  variables     T. 
with 


E[T  ]  ■     ,     "')    .     VarIT1]--Tra , 

J       p(yi  J  ,j  J       Pz(yrj+ir 


y-i-j+i 

1    -   -^7—  P 
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so  that 


and 


m  k"1      i 


p  >0  M  ' 


k-1 y-,    +  j 

1      r      P_ 


Var[Wk]  =  £    I    * 4 

k        P  j-0     (yrj)2 
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4.  Random  Search-Variant  2 

This  model  differs  from  Variant  1  only  in  that  the  cell  from  which  a 
ball  has  been  removed  is  not  selected  again.  The  development  parallels 
that  of  Variant  1  ^ery   closely  and  the  same  notation  and  definitions  as  in 
Section  3  are  used. 

The  basic  recurrence  for  the  random  variables  K(n)  ,  n  =  0,1,..., 
is  now 


K(n+1)  = 


Y]    -   K(n) 
K(n)     with  probability  1  -  -  _  K^  p  , 


y,  -  K(n) 
K(n)  +  1  with  probability  -  _  gr-x   p  , 


with  K(0)  =  0  .  To  avoid  the  trivial  case  yn  =  0  we  assume  that  yn  >  0  , 


0 


0 


i.e.     y,    <  m 


The  recurrence  for  the  distribution     Pn(k)     is  now 


Wk>  ■ 


yi  -  k 


l   -  — — r-  p  p  (k) 
m  -  k     H  rr  ' 


y     -  k  +  l 
+m-k+l     PpnCk-1)    . 


whence  the  generating  functions  G.  (t)  are 


sk(t) 


y(k) 

JkT(pt)k 


k 

n 

j-o 


1  -  t 


1  - 


y1  -  j 
m  -  J 


r  > 


U 
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and  calling 


yl  '  J 


p . 


we  have 


Thus 


y0p(j-i) 

aj  "  ai  '  (m-i)(m-j) 


A,  =  II  - 


y1  -  j 

m  -  j 


m  -  J 


m  -  i 


(ynP)  1-0 


k  *,n  J  -  i 


1  - 


yi   J  }     m(k)    M)k-j 


m  -  j 


(y0P 


k  ji(k-j)!  ' 


and  hence 


Gk(t) 


v=0 


>+k   [yl 


klk 
V 


I     MM- 


j=0 


1  - 


yl  -  J 

m  -  j 


v+k 


From  here  the  probabilities  are 


Pn(k) 


J_ 

k 
y0 


fy^  k 


j=o 


k-j 


1  - 


m  -  j 


-  i  ) 


k  =  0,...  ,y,  ,    n  =  0,1 , 
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Since  with  q  =  1  -  p 


1  - 


y-i  -  j 

m  -  J 


p  =  q  + 


m  -  j  ' 


we  can  also  write 


Pn(k) 


M     k 


3=0 


^  n 


q  + 


m  -  k  +  j 


k 


AkFn(0) 


by  using  the  identity 


AkF(x) 


I     (-l)k(k]F(x  +  k  -  j) 
j=0  lJj 


with 


F(x)   =  q  + 


m  -  x 


The  last  expression  can  be  used  to  obtain  the  factorial  moments  although 
the  results  are  somewhat  more  complex  that  for  Variant  1.  We  have  for 
r  -  0,...,y1 


y 


y(F)  =  I      k(r)P  (k) 


k=r 


y 


(r)  y-j-r 

I 


n+r 

y0  j=o 


y""j  Aj(ArFn(0)) 
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(r)  y.|-r 

I 


n+r 

y0  j=o 


j  j(y0-1)JArFn(n-j) 


by  using  the  identity 


b 

I 
j=0 


^)ab"jAjg(x)  =  I      (")(a-l)jg(x+b-j) 
JJ  j=o  UJ 


r,-n, 


with  b  =  y,  -  r  ,  a  =  yn  and  g(x)  =  A  F  (x)  . 


0 


Removing  the  difference  operator  and  substituting  for  F(x)  we  obtain  the 
formula 


(r)  . 


(r)  r  y-j-r 

I     I      M)1 


y;u  r  y] 


n+r 
yQ   i=0  j=0 


f*rrl         i 


YnP   ln 


0 


m+i+j-n-r 


r  =  0,...  ,y- 


In  particular,   the  expectation  is  (r  =  1) 


E[K(n) 


Vi     Jl 


yi-i 


yn      J=0 


fy-,-1 , 

V        J         J 


V 


q  +  — r*-i 

^       m-n+j-1 


q  + 


V 


m-n+j 


The  variance  Var[K(n)]  =  y^  +  E[K(n)]  -  E2[K(n)]  can  be  obtained 
similarly,  resulting  in  a  rather  long  expression. 
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Waiting  time  distribution  -  Variant  2 

The  derivation  is  analogous  to  that  of  Variant  1  with  the  only  difference 

yl  -  i  +  1 
being  the  parameter      H  +  ,  p  of  the  geometric  random  variables  T.  . 

The  probability  generating  function  <t>At)     of  the  waiting  time  W. 

is  now 


,(k) 


k-1 


>k(t)  =  -jjj  (Pt)  _nQ 


m 


1  -  t 


y-,   -  j 

1   - r-    p 

m  -  j  K 


-1-1 


so  that 


y1  -  k  + 


rXTT"PtGk-i^ 


y  tv+'<   pk 
v=0     ygdn-k+l) 


M   k-1  ,  ^id.jfk.! 


k 

i  j 


i  (-dk-'-j 


j=0 


y,  -  j 


>V 


+k-l 


P(Wk  -  n)  =  k 


_pk 


y0(m-k+l ) 


yl 

k 


k-1 

I  (-DJ 
j=0 


'k-1 


yop 


q  'h  m-k+l+j 


n-1 


n  =  k,  k+1 , 


k  =  1,...,  y. 
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From 


ftt   1   -  m.~  J   +  ] 
E[J]   "  P(yrj+D    ' 


Var[T.]   =-^1) 


yi  -  j  +  i 

m  -  j   +  1     K 


we  get  immediately 


and 


E[WR] 


,   k-1  . 

-  1    y    m  -  J    -  J< 

■  p  j£0  y]  -  J   "  p 


yo  k;]       l 


p  fiQ  y^  "  J 


k-l 


Var[Wk]   =  \    I 


m-j 


p'  j=q  yrJ 


i  - 


y^j 

m-j 


M    yo  k-i  (Hg)(yrj)  +y0 

2         2    .^n  ,        -x2 

P       p    j=o         (y^j) 
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5.  A  Multivariate  Search  Model 

In  this  section  we  investigate  some  aspects  of  the  multivariate  gener- 
alization of  the  Variant  1  model.  We  have  again  m  cells,  yQ  of  them 
empty,  and  m-yQ  of  them  containing  one  ball  each.  The  balls  can  be  of 
s  >_  1  different  types  and  initially,  we  have  y,  cells  containing  balls 
of  type  1,  y2  cells  containing  ball  type  2,  etc.  We  shall  refer  to  the 
empty  cells  as  containing  ball  type  0  ,  i.e.  no  ball.  Thus  the  initial 
configuration  is  specified  by  a  vector  of  nonnegative  integers 

y_  =  (y0>...,  ys)  , 

where  yQ  +  ...  +  y  =  m  and  y.  is  the  number  of  cells  containing  a  ball 
of  type  j  each. 

Cells  are  searched  in  repeated  trials  where  at  each  trial  a  cell  is 
selected  at  random  equally  likely  among  the  m  cells  and  independently 
of  previous  trials.  If  the  selected  cell  contains  a  ball  of  type  j  >  0 
the  ball  is  replaced  by  a  ball  of  type  j  -  1  .  Balls  type  0  are  just 
replaced,  i.e.  empty  cells  remain  empty.  The  same  cell  can  be  selected 
repeatedly  regardless  of  the  result  of  previous  trials.  Thus  for  s  =  1 
this  is  equivalent  to  the  Variant  1  model  with  p  =  1  .  The  model  is 
therefore  specified  by  a  single  vector  parameter  y_  . 
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Let  _X(n)  =  UQ(n),...,  X  (n))1  ,  n  =  0,1,...,  be  random  vector  where 


for  each  j  =  0, . . . ,  s 


th 


X.(n)  =  number  of  cells  of  type  j  after  the  n —  trial. 
Then  X(0)  =  y_  and  for  n  =  0,1,... 


Xn(n+1)  = 


X0(n) 


XQ(n)  +  1 


with  probability  1  - 


with  probability 


X-,(n) 


X^n) 


(5.1) 


Xj(n)  -  1 


X.(n+1)   U.(n)  +  1 


X-(n) 


with  probability 
with  probability 
with  probability  1 


m 
m 


Xj(n)  +  Xj+1(n) 


for  0  <  j  <  s  ,  and 


Xs(n+1)  =  { 


Xs(n)  -  1 


X  (n) 


with  probability 


with  probability  1  - 


X  (n) 

m   ' 


X  An) 


Of  course 


(5.2) 


XQ(n)  +  X-,(n)  +  ...  +  Xs(n)  =  m 
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for  all  n  =  0,1  ,. . .  . 
Alternately,  let 


K(n)  =  (KQ(n),...,  K  (n)) 


be  a  random  vector  where  for  j  =  0,...,  s  . 

K.(n)  =  number  of  trials  (up  to  n-th)  at  which  ball  type  j  was  found 

J 

Then  K(Q)  =  0  (zero  vector)  and 


f  yQ  +  K,(n) 

KQ(n)  +  1    with  probability     - 


KQ(n+l)  = 


m 


KQ(n) 


with  probability  1  - 


y0  +  ^(n) 


m 


y  +  K.+1(n)  -  K.(n) 
K.(n)  +  1    with  probability     -J JZJ J 


Kj(n+1)  =  { 


K,(n) 


with  probability  1  - 


y.  +  K.^Cn)  -  K.(n) 


for  0  <  j  <  s  ,  and 


K  (n)  +  1    with  probability     ys  "  Ks^n^ 


m 


Ks(n+1)  =  { 


Ks(n) 


with  probability  1  - 


yc  -  K  (n) 


s    s 


This  time  for  all  n  =  0,1,.. 


KQ(n)  +  ...  +  Ks(n)  =  n 
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There  is  a  simple  linear  relation  between  _X(n)  and  K_(n)  ,  namely 

XQ(n)  =  yQ  +  K^n)  , 
(5.3)         Xj(n)  =  yj  -  Kj(n)  +  KJ+1(n)  ,  0  <  j  <  s  , 

Xs(n)  =  ys  -  Ks(n)  , 
and  conversely  for  0  <  j   s 


Kj(n)  =  yj  +  ...  +  ys  -  (X..(n)  +  ...  +  Xfi(n))  , 


with 


K0(n)  =  n  -  (K^n)  +  ...  +  Ks(n))  . 

Note  that  in  terms  of  the  mine  sweeping  model  K, (n)  is  the  number  of 
deactivated  (exploded)  mines.  Also  B(n)  =  K,(n)  +  ...  +  K.(n)  ,  the 
total  number  of  balls  found,  is  the  total  number  of  contacts  the  sweeper 
makes  during  n  sweeps. 

Either  of  the  recurrence  relations  yields  immediately  the  corre- 
sponding recurrence  for  the  joint  probability  distribution  of  the  random 
vector  involved.  For  instance  denoting 

Pn(kr...,ks)  =  P(K1(n)  =  kr...,Ks(n)  =  k$) 

we  have  for  n  =  0,1 , . . .  . 
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Wl'—V  =^TTipn<lV---'ks' 


s-1  y,  +  k.+1  -  k.  +  1 


y,  -  k,  +  i 


s    s 


pn<kl ks-T  V'l 


with  initial  condition  PQ(k,,...,k  )  =1  if  k,  =  . . .  =  k  =  0  ,  and 
boundary  conditions  P  (k-,,...,k  )  =  0  whenever  k  <  0  for  some 

j  =  l,...,s  . 

In  principle,  this  can  be  used  to  recursively  evaluate  the  distribution 
P  (k-,,...,k  )  but  for  larger  values  of  s  this  is  not  practical  since  the 
storage  requirements  increase  exponentially  with  s  .  Generally,  an  array 
of  dimension  (y  +  1 )  (y  +  y     ^   +  1}  ...  (y  +  ...  +  y^  +  1)  would  be 
needed  to  store  the  current  values  of  P   although  some  savings  can  still 
be  made  due  to  the  fact  that  P  =  0  unless 

°ikj^j  +  kj+i  -  >iJ  <«  • 

We  therefore  attempt  first  to  evaluate  the  first  two  moments  of  the 
random  vector  _X(n)  . 

Define  the  multivariate  moment  generating  function 


*n(t)  =  E 


e*P  I     t  X  (n) 

1=0      l    l 


where  t_  =  (tQ,...,t  )  .  From  the  recurrence  (5.1) 
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exp  I     t  X  (n+1) 

1=0     l  l 


XQ(n)  +  I     X„(n)e 

U       £=0  l 


t   ,-t  ) 
2-1  I 


exp  I    t  X  (n) 


£=0 


£  £ 


Calling  temporarily 


1 

m 


if     £  =  0  , 


«t(t)  .  ^ 


l.**-!"**    if 

m 


0  <  £  <  s  , 


and  taking  expectation 


(5.4) 


*n+lte)  =  J0"t<D  af^nlD   • 


Now 


3Vh^         s     3a. (t)    3'*n^ 


s  324>   (t) 

+  Jo  ■»lt)  ^7 


whence  by  setting  t  =  0  . 
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S   3a. (0) 

E[X.(n+l)]  =  I    -r| ElX.(n)] 

1  £=1   3ti       l 


+  i  I  E[X.(n)X  (n)]  . 

m  1=0         1    * 


s  s 

However  I     E[X.(n)X£(n)]  =  E[X.(n)  J  X£(n)]  =  m  E[X.(n)] 


£=0 


z=0 


since   J  X  (n)  =  m  • 

£=0  * 


Further  for  0  <  I  <   s 


m 


9^(0) 
3t, 


~-     if     i  =  £  , 


<i 


if     i  =  £  -  1 


otherwise  , 


so  that 


E[XQ(n+l)]  =  E[XQ(n)3  +  jjEIX^n)] 


E[X.(n+l)]  =  (1  -  ^)E[X.(n)J  +  ^  E[X.+1(n)J 


for  0  <  i  <  s  ,  and 


E[Xs(n+l)]  ■  (1  -  ^)E[Xs(n)] 


Denoting  _X(n)  the  column  subvector 


(X1(n),...,Xs(n)) 
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and  by  Q  the  s  *  s  matrix 


Q  =  tqjJJ] 


with 


,P.) . 


1  -  - 
m 


m 


if     i  =  j  , 


if     j  =  i  +  1  , 


otherwise  , 


we  can  write 


whence 


E[X(n+l)]  =  QE[X(n)]  , 


E[X(n)]  =  Qn  y_ 


with  y_  =  (y-|,...,ys)  . 


However,  Q   is  an  upper  triangular  matrix  with  entries 


!  [  n  ' 
>i, 


m 


>)  = 


'ij 


(m-l)n-j+i     if     j   i 


if     j  <  i  , 


as  is  readily  verified  by  induction.  Hence 
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EUt(n)] 


1 

n  .L . 


J-i 


(m-1) 


n-j+i 


(5.6) 


1   p  fn)  /  i\n-£ 


i  =  l,...,s  ,  is  the  desired  expectation  vector, 


The  expectation  of  XQ(n)  is  obtained  from  (5.2) 


s  s-i 


^-^JhU")^11"^ 


s-l 


where  Z  =y  +  . . .  +  y  ,  a  *  1 , . . .  ,s  . 

The  expectation  of  the  random  vector  K(n)  is  obtained  using  the  linear 
relation  (5.3), 


EtV">j  =  v  3r  ?  SHn-<  v+1 


for     0  <  j   <  s   ,  and 


E[K°(n,1^-^i("Hn~'z^ 
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Recurrence  equations  for  the  second  moments  can  be  developed  similarly, 
Taking  second  partial  derivatives  of  (5.4)  we  get 


,K      3t-3t.  ~W 


3  X.  •  3  L  • 

1   J      Jl  =  l     1   J     I 


S   3a£(t)  3  -4>n(t) 

£=i   3ti    3tA 


S  3a   t   3  in  t 

+  y  l  ~ n  ~ 

A  3t.    3t-3t 

£=1  j    i  £ 


S         3\(t) 
^0  l~      3ti3tjH£ 


Again  upon  setting  t  =  0  the  last  terms  becomes  just  E[X.  (n)X.(n)  ]  and 
using  (5.5)  and 


(     1 
m 

if             i   = 

\iS)   _ 

"  m 

3t.3t.           ' 
1      J 

0 

otherwise   , 

i  =  j  =  i     or  i  =  j  =  i   -  1  , 


where  0<i<_j<_s  ,  0  <  I  <_  s  ,  we  obtain  for  0  <  i  <_  j  <_  s 
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ECX.fn+UXjfn+l)]  =  E£X1  (n)Xj(n)] 


-  -  E[X2(n)]   -  I  E[X2(n)] 

+  -  E[X.(n)X.^,(n)]  +  -  E[X . (n )X .    . (n )] 


^ECX^n)]  ♦iE[X1+l(n)]        if        j  -  1    , 


m 
-  -  E[X.(n)]    if    j  =  i  +  1  , 


if    j  >  i  +  1  , 


with  terms  E[X  (n)X  ,(n)]  and  E[X  ,(n)]  being  interpreted  as  zero. 
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This   again   looks   better  in  matrix  notation.      Denote 
Mn  -  [BlJ(n)]  ,  A-[0lJ]  ,  En  =[Eij(n)] 

the  s  x  s  matrices  defined  by 


Mi(j(n)  =  ECX^njX^n)]  , 


if     J  =  1  , 


aij  -\h      if      j  • '  + 1 


otherwise  , 


-  E[X.(n)]  +  i  E[X.  .(n)]     if 

^ECXj(n)]     if     j  -  1  +  1  , 
•-^EpC^n)]     if     1  =  j  +  1  , 


otherwise. 


Then  for  n  =  0,1,. . . 


J  =  i  , 


(5.7)  M  _,_.  =  M  +  AM  +  M  AT  +  E   , 

v    ;  n+1    n     n    n      n  ' 


with  T  denoting  transposition. 

Since  the  entries  of  Mn  are  just  y.y.  while  the  entries  of  E 

0      u     J  iJ  j  n 

are  available  from  (5.6)  for  all  n  =  0,1,...,  the  second  moment  matrix 

M„  can  be  evaluated  iteratively. 

n  J 


33 


It  is  even  possible  to  write  a  closed  form  expression  since  upon  iterating 
(5.7)  we  obtain 


(5.8) 


n   v 


Hn-  i   i  0(>\(A-y 

n   v=0  Z=0  v  ^    U 


k=0  v=0  £=0  V  *    n  i  K 


£     ( i) 
n  =  1,2,...  .  The  £-th  power  A  =  [a..]  of  the  matrix  A  is  an 

upper  triangular  matrix  with  entries 


(-1 


J"! 


J-i 


if 


J  >  1  , 


if 


j  <  i  , 


so  that  all  the  terms  on  the  right-hand  side  of  (5.3)  are  available. 

Having  the  mean  vector  and  the  second  moment  matrix  one  can  calculate 
the  covariance  matrix,  from  which  the  covariance  matrix  of  the  random 
vector  J<(n)  is  obtain  by  using  the  linear  transformation  (5.3).  Although 
having  an  expression  for  the  mean  vector  and  the  covariance  matrix  is  use- 
ful one  would  like  to  obtain  a  formula  for  the  probability  distribution  as 
we  did  in  the  univariate  case.  Unfortunately,  the  multivariate  case  is 
considerably  more  complicated  and  only  a  partial  result  is  obtained. 
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Let  0  <  j  <_  s  be  fixed  and  denote  by  Xv^(n)  the  number  of  cells 
among  the  ori  ginal  y.  eel  Is  (with  ball  type  i  initially)  that  after  n 
trials  contain  balls  of  type  j  .  Clearly 


X.(n)  =  x(0)(n)  + 
J      J 


.  *xjs>(n)  , 


where  of  course  X.  (n)  =  0  for  i  <  j  . 


Let  us  now  condition  on  the  random  vector  N  =  (N  ,...N  )  ,  where  N. 

—    o    s  i 

is  the  number  of  trials  which  resulted  is  a  selection  of  a  cell  from  the 
original  group  of  y.  cells.  Then  _N  is  multinomial  with  parameters  n 
and  y./m  ,  i  =  0,...,s  ,  and  X^.  '(n)  ,...,  X^.s'(n)  are  conditionally 
independent  given  N  .  It  follows  that 


(5.9) 


n     n 
o     s 

„    \  7o  '"h 
J         To--'nsJ    mn 


n 


P(X(.0)(n)  =  k|N)  *  P(X(.1)(n)  =  k|N)  *. 

J  J 


*  P(X<:s)(n)  =  k|N)  , 


where  the  summation  is  over  all  nonnegative  integers  n0>'-*>ns  sucn  that 

n  +  ...  +  n  =  n  and  asterisks  denote  convolutions. 
o        s 


Now  for  i  <  j  we  have   trivially 


P(X(.1  ^(n)  =  0|N)  =  1 

J 
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whi le  for  j  <  i  <  s  the  conditional  distribution  P(X .  (n)  =  k | H  =  n) 
is  identical  with  the  classical  occupancy  distribution  P^'(k)  of  findiny 


exactly  k  cell  occupied  with  exactly  i  -  j  balls  after  distributing  n- 
balls  randomly  into  y  -  y.  initially  empty  cells. 

This  distribution  is  most  easily  derived  from  the  generating  function 


CO       00 


H<y>(Uiv)  -  i    i  i£fi-vkP<*>(k) 


n=0  k=0 


by  conditioning  on  the  number  of  balls  placed  in  each  cell.   (See  [JK], 
p.  116  for  details.)  It  follows  that 


Hjy)(u'v)  =  M1^'^ 


and  since     P^    '(1)   =1     if     n  =  i -j     while     P^    '(0)   =  1     otherwise  we  have 
n  n 


5.10' 


H^y)(u,v)    - 


(v-1) 


,1  -J 


iy 


TT=J]T 


+  e 
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Now  for  0  <  k  <  y 


3kH^} 


dv 


(k)/  ui_j 


v=0 


i-J 


e  - 


mm 


y-k 


y-k 


.T-J 


k+i 


■'(k,X<-»T%hr'    * 


(y-k-i)u 


y-k 


=  y(k)  I  I  (-i)TO   (/-^)    d-J)(kn)tvt 

v=0  £=U     v  *  7  v!((1-j)!)K  * 


so  that  P^'(k)  -  — -  x  coefficient  at  u 
n        n 


Thus  with  y  =  y.  and  n  =  n.  we  have 


p.     (k)  =-J-rr-L-  i    (-D  (  \  ; 


ni 


n  . 


i   2=0 


(y-k+& 


ni  -  (i-j)(k+A) 


((i-j)!)k+£(n.  -  (1-j)(k+i) 


Upon  substitution  into  (S.9)  we  obtain  the  desired  expression  for  the 
probabilities  (for  0  <  j  :  s) 
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(y0  +  ...  +y,  ,)n  -  (nJ  +  -  +ns 

P(x.(n)  -lc)  =^-l  I  — 2 , JpL m 

J  mn  (n   -   (n^   +   ...    +  n$))! 


\  '       (-D  1  (£rVV 

B  ^n     r  !  (y  .-k.-i. ) !  k.+A. 

VU       '        '      '      '        ((i-j)!)   '      1(nr(i-j)(ki+4i))! 


where  the  first  two  sums  are,  respectively,  over  all  nonnegative  integers 
n.,...,n   and  k,,...,k   such  that 

n .+...+  n  <  n  , 
J         s  - 


and 


k .+...+  k  =  k  . 


Unfortunately,   for   larger  values   of     s,   k,   or  n     the  number  of  terms 
involved   in  the  indicated  sums   is   rather  prohibitive.     The  problem  of 
developing  a  manageable  computational   algorithm  to  evaluate  the  probabili 
ties     P(X.(n)   =  k)     is   left  for  future   investigation. 

It  should  be  pointed  out,   however,   that   the  generating  function 

oo  oo  n 


G,(u,v)   =     I       I    li™LvkP(X   (n)   =k) 
J  n=0  k=0  J 


is   just  a  product  of  the  generating  functions    (5.10),    specifically 


(yn+...+y  .    -,  )u      (v.)  (y    ) 

Sj(u.v)   =  e     U  J    l       H .   J    (u,v)...H.   s    (u,v)   . 
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Thus,  if  an  algorithm  can  be  found  to  expand  the  right-hand  side  of  this 
expression  into  a  power  series  in  u  and  v  the  probabilities  could  be 
determined  directly  from  this  expansion. 

As  an  example  of  this  approach  consider  the  case  when  the  parameters 
y,,...,y       are  themselves  random,  having  a  multinomial  distribution  with 
parameters  s  ,  b  =  y-|+...+y  ,and  p.  =  1/s  ,  i  =  l,...,s  ,  respectively. 
In  other  words  assume  that  each  ball  type  0  <  j  <_  s  was  selected  independ- 
ently and  at  random  among  all  nonzero  ball  types  l,...,s  .  Their  using  the 
multinomial  theorem  we  see  that 


G^u.v)  =  e 


V 


beu+(v-l)  \j^ 


ib 


Then  taking  k-th  partial   derivative  with   respect   to     v     and  setting 
v   =  0     we  get 


b-k  0/s-l     i 

Z=Q  \i=0 


k+£ 


.  (b-k-A)     (m-k-JOu 
b  e 


After  expanding  the  exponential    function   it   is   seen  that   the  coeffi 
cients   of  powers   of     u     will    be  expressible  as   ordinary  double  sums, 
involving  coefficients   of  powers   of  the  polynomial 


s-1     i 

y  — 

i  =0  ]  • 


The  latter  could  possibly  be  precomputed  and  thus  there  is  a  hope  for 
a  reasonably  efficient  algorithm. 
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